undef("calculate_mci")
function calculate_mci(inpath, outpath, stations, lats, year, MonthStart, MonthEnd);stations是站号,lats是纬度,year是年(注意应该大于等于1987),monthstart是起始月份(即每月最后一天),monthend是结束月份

;功能:计算气象综合干旱指数MCI(基于国标和指定日期过去30年的SPI\WPI)
;作者:王荣江
;相关论文:王荣江,李谢辉,周任君,王磊.三种气象干旱指数在四川省的适用性分析[J].长江流域资源与环境,2021,30(03):734-744.
;相关国标文件:《GBT 20481—2017  气象干旱等级》

;注意:处理含有缺测值\缺省值的数据时默认降雨量是0！气温为过去一周的平均温度！
;注意:站号是对应的气象站温度和降水文件，该文件格式默认为范例文件
;注意:这里需要调用该函数,在主文件main.ncl运行，做的循环比较多，处理时间可能会比较慢！见谅！

begin

    homedir = inpath ;这是数据读取路径
    mcis = new((/MonthEnd-MonthStart+1, dimsizes(stations)/), "float")

    do m=0, dimsizes(stations)-1

        data_rain = asciiread(homedir+stations(m)+"xxxx.txt", (/21915, 4/), "float"); 读取对应的站点气温数据,1957年-2016年所有日数
        data_temp = asciiread(homedir+stations(m)+"yyy.txt", (/21915, 4/), "float");读取对应的站点降水数据
        stanumber = stations(m); 站号
        nlat = lats(m);
        
        tt = 0;
        beginyear = 1987
        chooseyear = year
        beginyue = MonthStart
        endyue = MonthEnd
        lenyear = chooseyear - beginyear

        mci_1080 = new(endyue - beginyue + 1, float, -999);
        days = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
        ;if isleapyear(chooseyear) then
            ;days(1) = 29
        ;end if 
        if stanumber.eq.56152 then 
            days(6) = 30
            days(7) = 30
        end if 


        do i=0, 21914;
            if ismissing(data_rain(i, 3)) then 
                data_rain(i, 3) = 0.0
            end if 
            if (data_rain(i, 3).eq.-999) then
                data_rain(i, 3) = 0.0;   对缺测值\缺省值进行雨量为0的处理
            end if
        end do 

        do i=7, 21914;
            if ismissing(data_temp(i, 3)) then 
                data_temp(i, 3) = 0.0
            end if 
            if (data_temp(i, 3).eq.-999) then
                data_temp(i, 3) = avg(data_temp(i-7:i-1, 3)); 对缺测值\缺省值进行过去一周的平均温度处理
            end if
        end do 

        do years=lenyear, lenyear;
            do yue=beginyue, endyue;
                if yue.gt.12 then 
                    f = 12
                    years = lenyear + 1
                else
                    f = 0
                end if 
                do ri=days(yue-1-f), days(yue-1-f);
                    ;beginyear = 1987;
                    ;yue = 1;
                    ;ri = 1;
                    beginyear = 1987;
                    beginyear = beginyear + years;
                    a = 0.85;
                    x = 0.0;
                    print("the year is"+" "+beginyear);
                    print("the station is "+stanumber)
                    rain150 = new(30, float, -999);
                    rain90 = new(30, float, -999);
                    pre60 = new(61, float, -999);
                    rain30 = new(30, float, -999);
                    tempx = 0.0;
                    WAP = new(30, float, -999);

                    rainxx = new((/30, 3/), float, -999); 用于存储计算出的3个30年雨量序列
                    spi30 = new((/30, 3/), float, -999); 用于存储计算出的3个30年SPI化序列

                    do i=0, 21914;
                        if ((data_temp(i, 0).eq.beginyear).and.(data_temp(i, 1).eq.(yue-f)).and.(data_temp(i, 2).eq.ri)) then
                            tempx = avg(data_temp(i-29:i, 3));  获得过去30天的平均温度
                            break;            
                        end if
                    end do 

                    do j=0, 29;
                        do i=0, 21914;
                            if ((data_rain(i, 0).eq.beginyear).and.(data_rain(i, 1).eq.(yue-f)).and.(data_rain(i, 2).eq.ri)) then
                                rain150(j) = sum(data_rain(i-149:i, 3));获得过去30年间该日期前150天的降水总和数据  SPI150
                                rain90(j) = sum(data_rain(i-89:i, 3));获得过去30年间该日期前90天的降水总和数据    SPI90
                                rain30(j) = sum(data_rain(i-29:i, 3));获得过去30年间该日期前30天的降水总和数据
                                pre60 = data_rain(i-60:i, 3);
                            end if
                        end do 
                        xx = 0.0;
                        do k=0, 60;
                            xx = xx + a ^ (k) * pre60(k);  获得过去30年间该日期前61天的加权累积降水量 SPIW60
                        end do 
                        WAP(j) = xx;
                        beginyear = beginyear - 1;
                    end do 

                       
                    temp = new(12, float);
                    temp(0) = tempx;
                    temp(1:11) = 0;

                    pet = thornthwaite(temp, nlat, False, 0); 获得过去30天的蒸发量
                    pet@units = "mm/month";
                    MI = 0.0;
                    if pet(0).eq.0 then ;Hargreaves and Samani法
                        pet(0) = 0.0023 * (1368 / 2.45) * sqrt((max(temp)-min(temp))) * 17.8
                    end if 
                    MI = (rain30(0) - pet(0)) / pet(0);

                    rainxx(:, 0) = rain150; SPI150
                    rainxx(:, 1) = rain90;  SPI90
                    rainxx(:, 2) = WAP;    SPIW60

                    do k= 0, 2;
                        nyear = 30;
                        ;求尺度和形状参数
                            a1 = 0.0;
                            a2 = 0.0;
                            a3 = 0.0;
                            n1 = 0.0;
                            do i=0, nyear-1;
                                a1 = a1 + rainxx(i, k);
                                ;print(rainxx)
                                if ismissing(rainxx(i, k)) then 
                                    print(rainxx(i, k))
                                end if 
                                if (rainxx(i, k).gt.0) then 
                                    a2 = a2 + log10(rainxx(i, k))
                                else 
                                    n1 = n1 + 1;
                                end if 
                            end do
                            A = log10(a1 / (nyear - n1)) - a2 / (nyear - n1);
                            r = (1 + sqrt(1 + 4 * A / 3.0)) / (4 * A);
                            b = a1 / (nyear - n1) / r; 
                            ;调用函数计算GAMMA累积概率
                            p = new(nyear, float);
                            do i=0, nyear-1;
                                if (rainxx(i, k).gt.0) then 
                                    p(i) = rainxx(i, k) / b;高数换元法
                                    p(i) = gammainc(p(i), r);
                                else 
                                    p(i) = n1 / nyear;
                                end if 
                            end do 
                            ;对GAMMA分布概率进行正态化处理
                            c0 = 2.515517;
                            c1 = 0.802853;
                            c2 = 0.010328;
                            d1 = 1.432788;
                            d2 = 0.189269;
                            d3 = 0.001308;
                            do i=0, nyear-1;
                                if ((p(i).gt.0.5).and.(p(i).ne.1)) then 
                                    S = 1.0;
                                    p(i) = 1 - p(i);
                                    t = sqrt(log(1.0 / (p(i) * p(i))));
                                    spi30(i, k) = S * (t - ((c2 * t + c1) * t + c0) / (((d3 * t + d2) * t + d1) * t + 1.0));
                                else
                                    S = -1.0;
                                    t = sqrt(log(1.0 / (p(i) * p(i))));
                                    spi30(i, k) = S * (t - ((c2 * t + c1) * t + c0) / (((d3 * t + d2) * t + d1) * t + 1.0));
                                end if 
                            end do 
                    end do
                    a = 0.3;
                    b = 0.5;
                    c = 0.3;
                    d = 0.2;

                    MCI = 0.0;
                    Ka = 0.0;

                    if ((yue-f).eq.1) then
                         Ka = 1.0;
                    else if((yue-f).eq.2) then 
                        Ka = 1.0;
                    else if((yue-f).eq.3) then 
                        Ka = 1.1;
                    else if((yue-f).eq.4) then 
                        Ka = 1.2;
                    else if((yue-f).eq.5) then 
                        Ka = 1.0;
                    else if((yue-f).eq.6) then 
                        Ka = 1.2;
                    else if((yue-f).eq.7) then 
                        Ka = 1.2;
                    else if((yue-f).eq.8) then 
                        Ka = 1.2;
                    else if((yue-f).eq.9) then 
                        Ka = 1.2;
                    else if((yue-f).eq.10) then 
                        Ka = 1.0;
                    else if((yue-f).eq.11) then 
                        Ka = 1.0;
                    else if((yue-f).eq.12) then 
                        Ka = 1.0;
                    else 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    end if 
                    MCI = Ka * (a * spi30(0,2) + b * MI + c * spi30(0,1) + d * spi30(0,0));
                    mci_1080(tt) = MCI;
                    tt = tt + 1;
                end do 
                delete(rain150);
                delete(rain90);
                delete(pre60);
                delete(rain30);
                delete(WAP);
                delete(rainxx);
                delete(spi30);
            end do 
        end do 
        
        mcis(:, m) = mci_1080
        delete(mci_1080)
    end do 
    
    asciiwrite(outpath+stanumber+"_"+chooseyear+"_mci_test.txt", mcis);
    return(mcis)
end 